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ABSTRACT 

Radar  imaging  of  rotating  blade-like  objects,  such  as  helicopter  rotors,  using  narrowband 
radar  has  lately  been  of  significant  interest;  these  objects  cannot  be  adequately  described 
by  the  classic  point-scatterer  model.  Recently,  a  novel  ‘tilted-wire’  scatterer  model  has  been 
developed  that  can  provide  an  accurate  and  sparse  representation  of  radar  returns  from  such 
objects. 

Following  a  literature  review  on  compressed  sensing  algorithms,  covering  both  greedy  and 
lp  minimisation  methods  (0  <  p  <  1),  the  report  focuses  on  a  comparative  study  of  various 
greedy  pursuit  algorithms,  using  both  simulated  and  real  radar  data,  with  a  particular  em¬ 
phasis  on  the  use  of  the  tilted-wire  scatterer  model.  It  is  observed  that  the  greedy  algorithms 
that  select  multiple  atoms  at  the  matched-filtering  stage  do  not  perform  well  when  the  atoms 
used  in  the  dictionary  are  significantly  correlated.  Amongst  the  greedy  algorithms,  Ortho¬ 
gonal  Matching  Pursuit  (OMP)  exhibits  the  best  performance,  closely  followed  by  Conjugate 
Gradient  Pursuit  (CGP),  which  has  a  much  smaller  computational  complexity  than  OMP. 
In  applications  where  the  tilted-wire  model  requires  large  dictionaries  and  large  CPI  atoms, 
CGP  is  the  preferred  option. 
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A  Review  of  Sparsity-Based  Methods  for  Analysing 
Radar  Returns  from  Helicopter  Rotor  Blades 

Executive  Summary 

Signal  analysis  and  radar  imaging  of  fast-rotating  objects  such  as  helicopter  rotor  blades 
are  of  particular  research  interest  because  these  micro-Doppler  signals  cannot  be  processed 
by  conventional  range-Doppler  techniques  and  should  be  separated  from  other  non-rotating 
scattering  components  prior  to  further  processing.  In  addition,  the  point-scatterer  model, 
which  is  commonly  assumed  in  the  development  of  various  inverse  SAR  (ISAR)  and  radar 
tomography  approaches,  is  not  always  the  most  fitting  model  for  the  analysis  of  this  type  of 
micro-Doppler  signals.  The  novel  tilted-wire  model  offers  new  possibilities  to  overcome  the 
limitations  of  the  point-scatterer  model;  it  however  also  introduces  a  new  degree  of  complexity 
which  requires  the  use  of  state-of-the-art  sparsity-based  techniques. 

Since  the  tilted-wire  scatterer  model  can  be  used  to  facilitate  an  accurate  sparse  representation 
of  signals  from  rotating  blades,  a  comprehensive  review  of  known  algorithms  for  sparse  para¬ 
meter  estimation  techniques  is  carried  out,  covering  both  greedy  and  lp  minimisation  methods 
(0  <  p  <  1).  The  report  focuses  on  a  comparative  study  of  various  greedy  pursuit  algorithms 
using  both  simulated  and  real  helicopter  radar  data,  with  an  aim  to  accurately  estimate  the 
tilted-wire  parameters  associated  with  a  rotor  blade.  These  parameters  are  presented  as  scat¬ 
ter  plots  which  show  the  orientation,  length  and  tilt  of  the  estimated  wires  used  to  represent 
the  rotor  blade. 

Amongst  the  greedy  algorithms,  the  so-called  Orthogonal  Matching  Pursuit  (OMP)  technique 
exhibits  the  best  performance,  closely  followed  by  Conjugate  Gradient  Pursuit  (CGP),  which 
has  a  much  smaller  computational  complexity  than  OMP.  Important  improvements  and  ex¬ 
ploitation  of  these  modern  techniques  will  be  published  separately  in  the  near  future. 
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1  Introduction 


Compressed  sensing  (CS)  has  recently  received  remarkable  interest  thanks  to  its  wide  range  of 
applications  in  different  research  fields  including  electrical  engineering,  applied  mathematics 
and  computer  science.  In  radar  imaging,  CS  techniques  have  been  widely  exploited  in  many 
different  problems  such  as  sparse  phase  coherent  imaging  [1—3],  wide-angle  synthetic  aperture 
radar  (SAR)  imaging  for  anisotropic  scattering  [4-7],  multichannel  SAR  imaging  [8,  9],  and 
moving  target  indication  [10,  11].  In  this  report,  we  are  especially  interested  in  the  problem 
of  CS-based  radar  imaging  of  rotating  blade-like  objects  and  the  performance  of  different 
techniques  for  this  problem. 

Radar  signals  returned  from  rotating  blade-like  objects  are  usually  grouped  under  the  general 
category  of  ‘micro-Doppler’  signals,  where  the  Doppler  frequency  modulation  due  to  the  blade 
rotation  can  be  very  complex  and  extensive  [12-14].  Such  micro-Doppler  signals  can  not  be 
processed  by  conventional  range-Doppler  techniques  and  should  be  separated  from  other  non¬ 
rotating  scattering  components  prior  to  further  processing.  In  addition,  the  point-scatterer 
model,  which  is  commonly  assumed  in  the  development  of  different  inverse  SAR  (ISAR)  and 
radar  tomography  approaches,  are  not  always  the  most  fitting  model  for  the  analysis  of  micro- 
Doppler  signals. 

In  contrast  to  point-scatterers  with  isotropic  scattering  characteristics,  rotating  blade-like 
objects  behave  similarly  to  radial  antennas  with  a  lobe  structure  in  their  reflectivity  pattern. 
All  conventional  radar  imaging  theories  based  on  the  ideal  point-scatterer  are  thus  not  directly 
applicable.  A  plausible  extension  to  the  point-scatterer  model  is  the  ‘tilted-wire  scatterer’ 
model,  where  the  basic  scatterer  is  modeled  as  a  uniform  straight  wire  characterised  by  the 
position  of  its  centre  and  its  shape  [13,  14].  The  centre  position  vector  may  be  in  2D  or  3D 
space,  while  the  shape  parameters  include  at  least  the  finite  length  and  the  tilt  angle  of  the 
wire  relative  to  the  radial  direction. 

Rotating  blades,  in  particular,  can  be  approximated  as  a  collection  of  such  scatterers,  motiv¬ 
ated  by  the  fact  that  the  blades  may  exhibit  complex  bending  and  twisting  during  high-speed 
rotational  motion,  and  thus  a  strictly  radial  wire  model  may  not  be  sufficient  to  character¬ 
ise  the  scattering  from  the  blades.  More  importantly,  the  tilted-wire  model  facilitates  the 
application  of  CS  techniques  to  the  problem  of  radar  imaging  of  rotating  blades  [14]. 

The  fundamental  problem  in  compressed  sensing  is  to  recover  an  unknown  vector  x  E 
from  a  small  number  of  noisy  linear  measurements 

y  =  [y1,...,yM]T  =  ®x  +  n  E  CM,  (M  <  N).  (1) 

Here,  <l>  E  <CMxN  is  a  known  sensing  matrix  with  columns  normalised  to  unity  and  n  is  a  noise 
vector  with  energy  bound  of  ||n||2  <  e.  Since  the  sensing  matrix  $  provides  an  over  complete 
basis,  a  unique  solution  cannot  be  determined  using  the  conventional  inverse  transform  of  <I>. 
However,  if  x  is  sparse  or  compressible  (i.e. ,  well-approximated  as  being  sparse),  CS  theory 
enables  the  recovery  of  x  from  very  few  measurements  in  an  effective  and  robust  manner  [15— 
24].  The  sparse  recovery  problem  is  formulated  as  finding  the  sparsest  solution  of  x: 

x  =  argmin  ||cc||o  subject  to  x  E  B(y)  (2) 

xECn 
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where  the  constraint  x  €  B(y)  ensures  that  the  solution  x  is  consistent  with  the  measure¬ 
ments  y.  If  the  measurements  are  noise-free,  we  can  set 

B(y)  =  {x  :  $x  =  y}.  (3) 

Otherwise,  if  the  measurements  are  contaminated  by  bounded  noise,  we  can  instead  set 

B(y)  =  {x  :  \\y  -  <f>x\\2  <  e}.  (4) 

Solving  the  lo  minimisation  problem  (2)  is  NP-hard,  i.e.,  a  highly  non-convex  combinatorial 
optimisation  with  exponential  complexity,  and  thus  is  computationally  intractable  for  practical 
applications  [15-24].  An  attractive  alternative  is  to  consider  an  l\ -minimisation  problem 

x  =  argmin  ||a;||i  subject  to  x£B(y),  (5) 

x&CN 

which,  in  contrast  to  (2),  is  a  convex  optimisation  problem  (given  that  B(y)  is  convex)  and  can 
be  solved  effectively  in  polynomial  time  with  standard  convex  optimisation  techniques  [20,  21, 
23].  Importantly,  it  is  well-known  in  the  literature  that  the  l\ -minimisation  (5)  is  equivalent 
to  the  io-minimisation  (2)  for  the  sensing  matrices  satisfying  the  so-called  restricted  isometry 
property  (RIP)  with  a  constant  parameter  [16-19,  23,  24].  A  matrix  $  satisfies  the  RIP  of 
order  K  if  there  exists  a  5k  &  (0, 1)  such  that 

(1  -  5a')  llxlli  <  Il$xll2  <  (!  +  &K )  1 1 x  1 1 2  (6) 

holds  for  all  A'-sparse  vectors  x,  i.e.,  vectors  with  at  most  K  nonzero  components. 

Although  the  l\ -minimisation  formulation  based  on  convex  optimisation  provides  a  power¬ 
ful  framework  for  developing  computationally  tractable  CS  algorithms,  the  complexity  of 
i  i -minimisation  is  still  prohibitive  for  real-time  applications  [24-26].  Iterative  greedy  al¬ 
gorithms  [23-32]  have  recently  received  significant  attention  as  attractive  alternatives  to  con¬ 
vex  optimisation  techniques  thanks  to  their  low  complexity  and  simple  geometric  interpreta¬ 
tion. 

The  main  principle  behind  the  greedy  algorithms  is  that  they  identify  the  support  of  x  iter¬ 
atively,  where  one  or  more  elements  are  selected  at  each  iteration  based  on  some  greedy  rules 
and  their  contribution  is  subtracted  from  the  measurement  y.  The  main  advantage  of  the 
greedy  algorithms  is  that  they  are  more  computationally  efficient  than  the  l\ -minimisation 
algorithms.  In  terms  of  performance,  some  of  the  greedy  algorithms  have  been  shown  to 
have  theoretical  performance  guarantees  that  are  comparable  to  those  guarantees  derived  for 
convex  i \  -norm  optimisation  approaches  [21], 

Although  the  CS  literature  has  primarily  focused  on  the  ij -minimisation  and  greedy  pursuit 
algorithms,  there  exists  another  class  of  sparse  recovery  algorithms  based  on  lp- minimisation 
with  p  <  1.  The  ip-minimisation  with  p  <  1  was  shown  empirically  in  [33]  to  produce  exact 
recovery  with  fewer  measurements  than  with  p  =  1.  It  was  also  theoretically  proven  in  terms 
of  the  RIP  of  $  that  the  sufficient  condition  for  exact  recovery  of  x  via  the  ip-minimisation  is 
weaker  for  smaller  p  [34],  Although  the  ip-minimisation  with  p  <  1  is  a  non-convex  problem 
and  thus  computationally  intractable  as  described  in  the  literature  [34],  the  ip-minimisation 
with  p  <  1  can  be  effectively  solved  via  iteratively  re-weighted  least  squares  (IRLS)  [34, 
35], 
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In  this  report,  we  first  review  these  three  classes  of  sparse  recovery  algorithms  based  on  greedy 
pursuit,  /i -minimisation,  and  Zp-minimisation  with  p  <  1.  We  then  focus  our  attention  on  the 
greedy  pursuit  algorithms  and  present  a  comparative  performance  study  of  these  algorithms 
in  the  particular  problem  of  analysing  backscatter  signals  from  rotating  blades.  The  report  is 
organised  as  follows.  A  survey  on  the  CS  literature  is  presented  in  Section  2.  The  comparative 
performance  study  among  various  greedy  pursuit  algorithms  for  radar  imaging  of  rotating 
blades  is  presented  in  Section  3.  The  report  concludes  in  Section  4  with  a  brief  summary  of 
the  main  findings. 


2  Summary  of  Compressed  Sensing  Algorithms 

This  section  summarize  all  three  major  classes  of  algorithms  known  as  greedy,  Zi-minimisation, 
and  Zp-minimisation  algorithms.  However,  only  the  greedy  algorithms  will  be  studied  and 
discussed  in  great  detail. 


2.1  Greedy  Algorithms 


The  greedy  algorithms,  which  are  conceptually  simple  and  fairly  straightforward,  identify  the 
support  of  the  unknown  sparse  vector  x  progressively.  These  algorithms  usually  start  with  an 
initial  estimate  of 

k[0]  =  0,  (7) 

i.e.  the  support  set  -  the  indices  of  nonzero  elements  -  of  the  initial  estimate  is  A  =  0, 
and  an  initial  residual  of 

r[°l  =  y  -  $  x[0]  =  y.  (8) 

At  each  iteration,  one  or  more  columns  of  $  are  selected  based  on  the  correlation  values 
between  the  columns  of  $  and  the  residual  r.  The  indices  of  the  selected  columns  are  then 
added  to  the  support  set  A,  and  the  estimate  x  as  well  as  the  residual  r  are  also  updated  ac¬ 
cordingly.  The  greedy  algorithms  repeat  this  procedure  until  a  stopping  criterion  is  triggered, 
discussed  in  more  detail  in  Section  2.1.7. 

A  number  of  variants  of  iterative  greedy  algorithms  have  been  published  in  the  literature 
including  the  basic  Matching  Pursuit  (MP)  [27],  Orthogonal  Matching  Pursuit  (OMP)  [28], 
Stagewise  OMP  (StOMP)  [25],  Stagewise  Weak  OMP  (SWOMP)  [29],  Generalised  OMP 
(gOMP)  [26],  Regularised  OMP  (ROMP)  [23,  31],  Compressive  Sampling  Matching  Pur¬ 
suit  (CoSaMP)  [32],  Subspace  Pursuit  (SP)  [24],  Gradient  Pursuit  (GP)  [30],  Conjugate  GP 
(CGP)  [30],  and  Stagewise  Weak  CGP  (SWCGP)  [29].  Some  of  these  variants  will  be  discussed 
in  further  detail  in  the  following  sub-sections. 


2.1.1  Basic  Matching  Pursuit  (MP) 

The  simplest  greedy  algorithm  is  MP  [27]  which  is  summarised  in  Algorithm  1.  At  each 
iteration,  MP  selects  one  column  of  which  yields  the  largest  correlation  with  the  current 
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Algorithm  1  MP  Algorithm 


1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 


procedure  INPUT:  y,  OUTPUT:  rW,  xU 
Initialisations:  r^  =  y,  x'°!  =  0. 

for  iteration  i  =  1;  i  :=  i+1  until  stopping  criterion  is  met  do 
Identify: 

C  —  ^7  T 

j*  =  argmaxj  |cj| 

Update: 

-W  -  [*— i]  , 

xu  =  xj  ■  +  cr 

rW  =  r[®_1]  _ 

end  for 

end  procedure 


residual,  and  updates  the  set  of  coefficients  x  accordingly.  Here,  the  coefficient  update 


(9) 


may  occur  for  a  coefficient  that  has  been  updated  in  an  earlier  iteration.  In  other  words, 
MP  may  select  a  particular  column  of  $  multiple  times  and  refine  its  coefficient  to  enhance 
the  approximation  performance.  The  computational  complexity  of  MP  in  each  iteration  is 
dominated  by  the  matrix  vector  product  in  the  identification  step  which  requires  a 

computation  of  0(MN )  for  unstructured  matrices  or  0(N  log(M))  for  structured  matrices 
by  exploiting  fast  Fourier  transform  (FFT)-based  algorithms  [21,  30].  In  addition,  finding  the 
largest  element  in  c  in  the  identification  step  and  updating  the  residual  rM  require  N  opera¬ 
tions  and  M  operations  respectively. 


2.1.2  Orthogonal  Matching  Pursuit  (OMP) 

The  OMP  algorithm  [28] ,  a  more  sophisticated  refinement  of  MP,  is  one  of  the  most  popular 
greedy  algorithms.  The  pseudo-code  of  OMP  is  summarised  in  Algorithm  2.  The  main 
difference  between  OMP  and  MP  is  that,  in  each  iteration,  OMP  projects  y  orthogonally  onto 
the  columns  of  $  associated  with  the  current  support  set  AW,  in  a  least-squares-error  sense, 
to  obtain  a  new  approximation  of  x.  In  other  words,  OMP  minimises  ||y  —  <f*xM  over  all  xM 
with  support  AW.  Another  difference  is  that  OMP  only  selects  an  element  maximally  once 
and  the  residual  r^  is  always  orthogonal  to  the  current  selected  element  set  3>A[i]  as  a  result 
of  the  least-squares  estimation.  This  least-squares  estimation  makes  OMP  computationally 
more  demanding  than  MP.  However,  OMP  provides  a  superior  performance  compared  to  MP, 
particularly  in  terms  of  its  convergence  property. 

The  computational  complexity  of  OMP  depends  on  the  actual  implementation  of  the  least- 
squares  estimation.  For  example,  the  QR  factorisation  approach  requires  2 Mk  +  3 M  opera¬ 
tions  to  obtain  the  least-squares  solution  where  k  is  the  size  of  the  current  support  set  [30]. 
It  also  depends  on  the  number  of  iterations  needed  for  signal  recovery,  i.e. ,  in  the  order  of 
O(KMN)  for  K  iterations  [24,  26,  31].  This  complexity  is  significantly  smaller  than  that 
of  the  l\ -minimisation  algorithms  based  on  convex  optimisation  [23,  24,  36].  However,  OMP 
does  not  offer  the  strong  theoretical  guarantees  as  the  L -minimisation  techniques  [23,  31,  36]. 
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Algorithm  2  OMP  Algorithm 


1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 


procedure  INPUT:  y,  OUTPUT:  rW,  xU 
Initialisations:  r^  =  y,  =  0,  =  0. 

for  iteration  i  =  1;  i  :=  i+1  until  stopping  criterion  is  met  do 
Identify: 

c  —  i] 


j*  =  argmaxj  |cj| 

Merge  Supports: 

AW  =  At*-1!  U  j* 

Update: 

=  argminx  ||y  —  x|| 2 ,  where  <&A[*]  are  at°ms  in  current  set  AW 

r^y-^AW*!^ 

end  for 

end  procedure 


It  was  demonstrated  theoretically  and  empirically  in  [37]  that  OMP  can  reliably  recover  a 
A-sparse  signal  for  given  M  =  0{K  log  N)  noise-free  measurements. 

Specifically,  for  a  Gaussian  sensing  matrix  (its  atoms  are  drawn  independently  from  the  stand¬ 
ard  Gaussian  distribution),  OMP  can  recover  any  A-sparse  signal  x  with  probability  exceeding 
(1  —  25)  given  that  5  6  (0,  0.36)  and  M  >  CKlog(N/5 )  [37].  The  constant  C  satisfies  C  <  20 
or  can  be  reduced  to  C  ~  4  for  a  large  value  of  A.  However,  in  contrast  to  ^-minimisation,  the 
recovery  guarantees  of  OMP  is  nonuniform  [23,  38,  39].  In  particular,  when  M  =  0(K  log  N), 
there  is  high  possibility  of  existence  of  a  Ji-sparse  vector  x  for  which  OMP  will  select  an 
incorrect  element  at  the  first  iteration  for  certain  random  matrices  <1?  [38,  39].  It  has  been 
shown  in  [38]  that  uniform  recovery  guarantees  for  OMP  are  impossible  for  the  natural  random 
sensing  matrix  3». 

In  addition,  the  conditions  on  the  sensing  matrices  required  by  OMP  are  more  restrictive  than 
the  RIP  condition  [23,  24],  In  particular,  Tropp  [40]  showed  that  the  sufficient  condition  for 
OMP  to  recover  a  A-sparse  x  exactly  from  noise-free  measurements  is 


A  < 


1 

2A  -  1 


(10) 


where  /u  is  the  mutual  coherence  parameter  defined  by  the  maximum  value  of  the  modulus  of 
the  inner  product  between  two  distinct  columns  of  the  sensing  matrix  $  (commonly  known  as 
the  Mutual  Incoherence  Property  (MIP)  introduced  in  [41]).  This  result  was  then  extended  to 
the  case  of  noisy  measurements  in  [42].  Note  that  if  MIP  holds  then  RIP  also  holds,  but  the 
converse  is  not  true.  Thus  the  MIP  condition  is  more  restrictive  than  the  RIP  condition  [42] . 
However,  the  advantage  of  the  MIP  over  the  RIP  is  that  the  MIP  is  more  straightforward  to 
verify  for  any  given  matrix  $  [42], 

Although  the  RIP  has  been  proved  in  [43]  to  hold  with  high  probability  for  the  random  con¬ 
struction  of  the  sensing  matrix  <J>,  there  is  still  probability  of  failure  of  the  RIP  for  a  particular 
realisation  of  a  random  matrix.  In  addition,  testing  whether  a  given  matrix  satisfies  the  RIP 
is  NP-hard  and  computationally  infeasible  [44].  In  terms  of  RIP,  Wakin  and  Davenport  [45] 
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Algorithm  3  StOMP  and  SWOMP  Algorithms 

1:  procedure  INPUT:  y,  OUTPUT:  rM,  xU 
2:  Initialisations:  r^  =  y,  x)°]  =  0,  =  0. 

3:  for  iteration  i  =  l;*:=*  +  l  until  stopping  criterion  is  met  do 

4:  Identify: 

5:  c  =  3>rr[*-1l 

6:  r  =  {j  :  \Cj\  >  AW} 

7:  Merge  Supports: 

8:  AW  =  A[*-11  U  J* 

9:  Update: 

10:  =  arg minx  ||y  -  $aWx||2 

11:  rW=y-^AWxWw 

12:  end  for 

13:  end  procedure 


showed  that  the  RIP  of  order  K  +  1  with 

*K+1  <  (11) 

is  sufficient  for  OMP  to  exactly  reconstruct  any  A'-sparse  signal  from  noise- free  measurements. 
This  condition  was  then  improved  to  [46] 

s,f+1 "  7W+i‘  (12) 

2.1.3  Stagewise  OMP  (StOMP),  Stagewise  Weak  OMP  (SWOMP)  &;  Gen¬ 
eralised  OMP  (gOMP) 

Although  OMP  is  computationally  cheaper  than  the  ^-minimisation  algorithms,  it  may  still 
not  be  the  best  choice  for  certain  class  of  problems,  especially  for  not  very  sparse  signals 
and  in  terms  of  computational  complexity  [25,  26,  29,  30].  The  main  reason  for  this  is  the 
fact  that  OMP  only  selects  a  single  column  of  the  sensing  matrix  $  (a  single  atom),  at  each 
iteration  and  thus  it  must  run  at  least  as  many  iterations  as  the  number  of  nonzero  elements 
in  the  solution.  This  computational  issue  motivates  the  selection  of  multiple  atoms  at  a  time, 
as  has  been  incorporated  in  a  number  of  variants  of  OMP  including  StOMP  [25],  SWOMP 
[29],  and  gOMP  [26].  Although  the  computational  complexity  for  one  iteration  of  StOMP, 
SWOMP  and  gOMP  is  similar  to  that  of  OMP,  the  StOMP,  SWOMP  and  gOMP  algorithms 
require  fewer  iterations  to  obtain  the  same  number  of  nonzero  elements  in  the  solution  when 
compared  to  OMP.  This  makes  StOMP,  SWOMP  and  gOMP  computationally  more  effective 
than  OMP. 

The  maximum  selection  operation 

AW  =  At*-1]  U  j* ,  with  j*  =  argmax|cj|  (13) 

j 

in  OMP  is  replaced  with  a  thresholding  operation  in  StOMP  and  SWOMP  as  follows: 

A[*]  =  A[i-1]  U  {j  :  |Cj-|  >  Al*]}  (14) 


6 


UNCLASSIFIED 


UNCLASSIFIED 


DST-Group-TR-3292 


Algorithm  4  gOMP  Algorithm 

1:  procedure  INPUT:  y,  L.  OUTPUT:  rW,  xU 
2:  Initialisations:  r^  =  y,  =  0,  =  0. 

3:  for  iteration  i  =  l;*:=*  +  l  until  stopping  criterion  is  met  do 

4:  Identify: 

5:  c  =  4>rr[*_1] 

6:  J*  =  {indices  of  the  L  largest  elements  of  c} 

7:  Merge  Supports: 

8:  AW  =  A[*_11  U  J* 

9:  Update: 

10:  x^w  =  arg minx  ||y  -  $aWx||2 

11:  r[,1=y-%li 

12:  end  for 

13:  end  procedure 


where  AW  is  the  threshold  at  iteration  i.  The  full  algorithms  of  StOMP  and  SWOMP  are 
summarised  in  Algorithm  3.  The  two  variants  differ  only  in  the  choice  for  this  threshold. 

The  threshold  used  in  StOMP  is  given  by 


' StOMP 


Vm 


(15) 


where  the  typical  value  of  the  parameter  fW  is  2  <  fW  <  3.  Two  threshold  strategies  of  f  W  were 
proposed  explicitly  in  [25]  for  the  case  of  random  sensing  matrix  $  generated  from  a  uniform 
spherical  ensemble,  i.e.  the  columns  of  $  are  independently  and  identically  distributed  (i.i.d.) 
points  on  the  unit  sphere.  These  strategies  were  motivated  by  classical  detection  criteria  of 
false  alarm  control  and  false  discovery  control.  However,  theoretical  performance  guarantees 
of  StOMP  are  not  available  for  more  general  sensing  matrices  <!>. 

Moreover,  from  the  practical  perspective,  the  choice  of  the  parameter  f M  appears  critical  to 
the  performance  of  StOMP,  but  the  question  on  the  optimal  selection  of  fM  has  not  been 
fully  addressed  in  the  literature  [29].  It  was  demonstrated  in  [29]  that  StOMP  sometimes 
terminates  prematurely  when  all  inner  products,  i.e.,  |cj|  for  all  values  of  j,  fall  below  the 
threshold  A^.  In  addition,  the  simulations  in  [29]  showed  mixed  results  on  the  performance 
of  StOMP. 

The  threshold  used  in  SWOMP  is  given  by 


^ SWOMP  ~  am4x{lcjl} 
j 


(16) 


where  a  6  (0, 1]  is  the  ‘weakness  parameter’.  It  is  important  to  note  that  SWOMP  enjoys  a 
weakened  version  of  the  performance  properties  of  OMP  [40]. 

In  contrast,  as  summarised  in  Algorithm  4,  gOMP  selects  a  fixed  number  of  atoms  for  each 
iteration,  i.e. 

AW  =  A[i-i]  u  {</*(!), . . . ,  J*{L )}  (17) 
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Algorithm  5  ROMP  Algorithm 


1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 
17 


procedure  INPUT:  y,  4>,  K.  OUTPUT:  rW,  xU 
Initialisations:  r^  =  y,  =  0,  Ai°]  =  0. 

for  iteration  i  =  1;  i  :=  i+1  until  stopping  criterion  is  met  do 
Identify: 

C  —  A 

J  =  {indices  of  the  K  largest  elements  of  c} 

Regularise: 

J*  =  arg  max  j;  ||cJ;  || 2 

where  Ji  is  all  subsets  of  J  with  comparable  magnitudes: 
|ca|  <  2|c;,|  for  all  a,  b  E  J; 

Merge  Supports: 

AW  =  At*-1!  U  J* 

Update: 

=  arg minx  ||y  —  ^A[i]x||2 

rt^y-^AwxJl, 

end  for 

end  procedure 


where  L  is  the  number  of  atoms  selected  at  each  iteration,  and  J*(l )  is  the  index  of  the  i-th 
largest  absolute  value  of  the  entries  of  c  explicitly  defined  as 

J*{1)  =  arg  max  |cj|.  (18) 

For  the  noise-free  measurement  scenario,  the  RIP  order  of  LK  with 

Slk<  Vk  +  3Vl’  {K>1)  (19) 

is  a  sufficient  condition  for  gOMP  to  obtain  an  exact  recovery  of  any  A-sparse  vector  within 
K  iterations  [26].  A  performance  bound  on  the  estimation  error  for  the  signal  reconstruction 
in  the  presence  of  noise  was  also  derived  in  [26]. 

2.1.4  Regularised  OMP  (ROMP) 

An  alternative  modification  on  the  identification  step  was  proposed  in  the  ROMP  algorithm  [23, 
31]  which,  at  each  iteration,  selects  K  largest  entries  of  the  correlation  c  and  groups  them 
into  subsets  J\  with  comparable  magnitudes: 

[ca|  <  2|cft|  for  all  a,  b  E  J/.  (20) 

ROMP  then  selects  the  set  Ji  with  the  maximum  energy  1 1  c  ^  1 1 2  •  The  algorithm  is  summarised 
in  Algorithm  5. 

In  contrast  to  OMP  and  its  aforementioned  variants,  ROMP  enjoys  uniform  recovery  guaran¬ 
tees  [23].  In  particular,  the  RIP  of  order  2 1\  with 

52K  <  0.03/Vlog  I<  (21) 
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is  the  sufficient  condition  for  ROMP  to  exactly  reconstruct  any  A"-sparse  signal  x  from  its 
noise- free  measurement  y  =  <l>x  [23].  For  the  scenario  with  the  presence  of  noise  n  in  the 
measurement  y  =  $x  +  n,  the  theoretical  performance  bound  of  ROMP  is  given  by 

||x  -  x||2  <  KM-y/logAT ||n||2,  (22) 

as  discussed  in  [31].  In  fact,  the  logarithmic  factor  in  this  expression  yields  a  stronger  the¬ 
oretical  requirement  for  the  RIP  compared  to  that  of  the  / 1 -minimisation  methods  [29,  31]. 
However,  from  a  practical  point  of  view,  one  may  be  more  interested  in  the  average  per¬ 
formance  than  the  worse  case  performance  provided  by  these  theoretical  analyses.  It  was 
empirically  demonstrated  in  [29]  that  the  average  performance  of  ROMP  was  notably  worse 
than  OMP  and  StOMP. 

The  computation  for  one  iteration  of  ROMP  is  slightly  more  demanding  than  that  of  OMP 
due  to  the  differences  in  the  identification  and  regularisation  steps.  The  selection  of  K  largest 
elements  of  c  can  be  done  via  a  sorting  algorithm  which  typically  requires  a  complexity  of 
O  (N  log  N)  [31].  Since  the  selected  support  set  J  is  already  sorted,  the  regularisation  step 
can  be  performed  effectively  by  searching  over  consecutive  intervals  of  J,  and  thus  it  only 
requires  the  complexity  of  0(K )  [31].  However,  ROMP  may  select  multiple  atoms  at  a  time. 
As  a  result,  ROMP  may  require  less  run  time  than  OMP  if  a  smaller  number  of  iterations  is 
required  for  ROMP  to  achieve  the  same  sparsity  level  (i.e.,  same  number  of  non-zero  elements) 
in  the  solution. 


2.1.5  Compressive  Sampling  MP  (CoSaMP)  &  Subspace  Pursuit  (SP) 


In  contrast  to  the  other  aforementioned  variants  of  OMP  which  only  focus  on  the  modifica¬ 
tion  to  the  identification  step  of  OMP,  the  CoSaMP  algorithm  proposed  in  [32]  and  the  SP 
algorithm  proposed  in  [24]  include  an  additional  pruning  step  in  each  iteration.  The  CoSaMP 
and  SP  algorithms  are  summarised  in  Algorithm  6.  The  main  idea  behind  CoSaMP  and  SP  is 
that  they  maintain  a  fixed  number  of  nonzero  elements  in  each  active  set  A^  for  each  iteration 
by  removing  insignificant  elements  via  the  pruning  step. 

Note  that,  in  the  other  greedy  algorithms  previously  described,  once  an  atom  is  selected,  it 
will  always  stay  in  the  active  set  until  the  algorithm  terminates.  Specifically,  the  pruning 
step  in  CoSaMP  and  SP  only  retains  K  largest  entries  in  the  least-squares  optimisation  for 
the  merged  support  set  A W .  CoSaMP  and  SP  then  re-estimate  the  least-squares  solution 
corresponding  to  the  retained  support  set  AW.  The  only  difference  between  CoSaMP  and  SP 
is  that  CoSaMP  adds  2 K  new  atoms  to  the  active  support  set  in  the  identification  step  while 
SP  only  adds  K  new  atoms. 

The  main  advantage  of  these  two  algorithms  is  that  they  are  not  computationally  complex 
and  provide  strong  theoretical  guarantees  that  are  comparable  to  those  derived  for  the  convex 
/i-optimisation  methods  [21,  24,  32],  However,  the  drawback  of  CoSaMP  and  SP  is  that  they 
require  the  knowledge  of  the  sparsity  level  K  as  an  input,  which  may  not  be  known  a  priori 
in  practice. 
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Algorithm  6  CoSaMP  and  SP  Algorithms 

1:  procedure  INPUT:  y,  K.  OUTPUT:  rW,  xU 
2:  Initialisations:  r^  =  y,  =  0,  =  0,  A^l  =  0. 

3:  for  iteration  i  =  \]  i  :=i  +  \  until  stopping  criterion  is  met  do 

4:  Identify: 

5:  c  =  $Tr[*_1l 

6:  J*  =  {indices  of  the  2 K  largest  elements  of  c}  (CoSaMP) 

7:  or  J*  =  {indices  of  the  K  largest  elements  of  c}  (SP) 

8:  Merge  Supports: 

9:  AW  =  A[*-11  U  J* 

10:  Estimate  Least  Squares  (LS)  Solution: 

11:  x^w  =  arg minx  ||y  -  $A[i]x||2 

12:  Pruning: 

13:  AW  =  {indices  of  the  K  largest  elements  of  xj^} 

14:  Update: 

15:  x^.,  =  arg minx  ||y  -  $A[i]x||2 

16:  r[i]  =y-^AWiAH 

17:  end  for 

18:  end  procedure 


2.1.6  Gradient  Pursuit  (GP),  Conjugate  GP  (CGP)  &  Stagewise  Weak 
CGP(SWCGP) 


Another  trend  of  greedy  pursuit  technique  based  on  directional  updates  was  proposed  in  [29, 
30],  namely  GP,  CGP  and  SWCGP.  The  main  idea  behind  these  directional  algorithms  is 
that  they  exploit  directional  optimisation  to  update  the  coefficients  of  the  selected  elements 
in  each  iteration  instead  of  using  orthogonal  projection,  i.e.,  least-squares  estimation,  as  in 
OMP  and  its  variations.  The  replacement  of  the  costly  orthogonal  projection  by  the  directional 
optimisation  leads  to  computationally  more  efficient  algorithms  while  still  retaining  similar 
performance  as  OMP.  Directional  optimisation  refers  to  an  iterative  technique  of  finding  a 
local  minimum  of  a  given  cost  function  by  starting  with  an  initial  point  in  the  parameter 
space  and  moving  towards  the  direction  that  minimises  the  cost  function.  A  great  example 
of  directional  optimisation  is  the  gradient  method  where  the  update  direction  is  determined 
by  the  gradient  of  the  cost  function  at  the  current  point. 

The  expression  of  directional  update  is  [29,  30] 


=  x 


[<-l] 

'AW 


+  <AJd 


am 


(23) 


where  is  the  update  direction  and  aN  is  the  step  size.  The  optimal  value  of  the  step 
size  to  minimise  ||y  —  <1>xM|||  over  all  xM  with  support  A^  (i.e.,  the  same  quadratic  cost 
as  in  OMP)  is  explicitly  given  by  [30,  47] 


(rH,  bW) 


(24) 
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Algorithm  7  Gradient  Pursuit  Algorithm 


1: 

2: 

3: 

4: 

5: 

6: 

7: 

8: 

9: 

10: 

11: 

12: 

13: 

14: 

15: 


procedure  INPUT:  y,  OUTPUT:  rW,  xU 
Initialisations:  r^  =  y,  =  0,  =  0. 

for  iteration  i  =  1;  i  :=  i+1  until  stopping  criterion  is  met  do 
Identify: 

C  —  ^7  1] 


j  =  argmaxj  |Cj| 

Merge  Supports: 

AW  =  U  j* 

Directional  Update: 

compute  d^]  (see  text  for  details) 
bW  =  $,.,id[*L  n.W  =  <rl'l;.blij) 


x 


AH 


■  Ai 

_ 

AW 


dAH>  flL‘J  = 


=  X 


+  a  I 


i]d[ 


l|bW||l 


ri-j  =  r [®— i]  _  aWb[®] 


AH 


end  for 

end  procedure 


where  bW  =  $A[i]d^[i]  and  (•)  denotes  the  inner  product  operation.  The  pseudocode  for 
directional  pursuit  family  is  given  in  Algorithm  7.  Note  that  MP  and  OMP  can  be  naturally 
cast  into  the  framework  of  directional  pursuit  with  update  directions  of  5^  and  <&ahcaW  (the 
superscript  f  denotes  the  matrix  pseudoinverse  operation). 

GP  utilises  the  negative  gradient  of  the  cost  quadratic  cost  function  ||y  —  <&xFl  |||  as  the  update 
direction  [30] 

dAH  =  *AW  (y  -  #AH  XAw1] )  =  CAW  (25) 

where  cA[ij  is  a  subvector  of  c  with  support  AW. 


On  the  other  hand,  the  directional  update  in  CGP  exploits  the  conjugate  gradient  method 
which  is  widely  used  to  solve  quadratic  optimisation  problems  [47].  Fundamentally,  to  min¬ 
imise  a  cost  function  of  (l/2)x2  Gx-yJx  (he.,  equivalent  to  solving  y  =  Gx  for  x),  the  con¬ 
jugate  gradient  method  successively  applies  line  minimisations  along  directions  which  are  G- 
conjugate,  where  set  of  directions  {dW .  dl2! , . . . ,  dM  }  is  defined  as  G-conjugate  if  (d^ ,  Gd^l)  = 
0  for  all  i  A  j.  Note  that  (•,  •)  denotes  the  inner  product  between  two  vectors.  In  the  frame¬ 
work  of  directional  pursuit,  CGP  aims  to  minimise  the  cost  function  of  ||y  —  <f,A[l]XA[l]  |2  by 
calculating  an  update  direction  that  is  GA[i]-conjugate  to  all  previous  update  directions.  Here, 

Gah  =  *ah  ■ 

The  main  advantage  of  the  conjugate  gradient  method  is  that  the  conjugate  directions  can  be 
computed  iteratively  using 

dW  =  cw  +  (26) 

if  the  first  conjugate  direction  is  initialised  to  dM  =  —  cM.  Here,  [3^  is  given  by  j3^  = 
(cW,  Gd[*_1l)/(d[*_h,  Gdl!_1l)  to  ensure  (dH,GdI®_1])  =  0.  It  is  important  to  note  that  this 
procedure  leads  to  d^  conjugate  to  all  previous  directions  dM, . . .  ,  d^-1!.  This  principle  is 
adopted  to  compute  the  update  direction  of  CGP,  i.e. , 


(27) 
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where 


(3® 


(  (&AH-1]  dJj[4S] ) ,  (*AW  caM  )  ^ 


«AMd 


I*-1]  112 

At*— i]  N  2 


(28) 


However,  in  contrast  to  that  in  (26),  the  iterative  procedure  in  (27)  only  guarantees  the 
conjugacy  of  to  the  last  preceding  direction  but  not  to  all  the  previous  directions 

because  the  dimensionality  of  the  solution  space  varies  over  iterations.  Therefore,  using  (27) 
leads  to  an  approximate  version  of  CGP  (referred  to  as  CGP  throughout  this  report  for 
simplicity).  The  more  sophisticated  version  of  CGP  can  be  found  in  [29],  where  each  new 
direction  is  calculated  to  ensure  the  conjugacy  with  all  previous  directions.  However,  this 
version  is  computationally  more  complex  than  the  approximated  CGP. 

In  common  with  MP,  the  identification  step  of  GP  requires  a  matrix  vector  multiplication 
and  N  operations.  Computing  the  size  step  for  the  directional  update  of  GP  requires  an 
extra  matrix  vector  multiplication  and  2 M  operations  [30].  In  addition,  updating  takes 
k  operations  where  k  is  the  size  of  the  current  support  set,  and  updating  the  residual  rW 
can  be  done  via  M  operations.  On  the  other  hand,  compared  to  GP,  CGP  requires  extra 
computations  to  compute  the  update  direction  including  one  additional  matrix  vector 
multiplication  and  M  +  k  additional  operations  [30]. 

Similar  to  MP  and  OMP,  GP  and  CGP  only  selects  a  single  element  at  each  iteration  making 
it  not  suitable  for  large-scale  problems  in  terms  of  computational  performance.  Motivated 
by  this  issue,  a  stagewise  weak  version  of  CGP  (SWCGP)  was  proposed  in  [30]  allowing 
multiple  elements  to  be  selected  at  each  iteration.  Since  the  extension  of  CGP  to  SWCGP  is 
analogous  with  the  extension  of  OMP  to  SWOMP  which  has  been  discussed  in  Section  2.1.3, 
the  details  of  SWCGP  is  omitted  here.  The  interested  readers  may  refer  to  [30]  for  more 
detailed  description  of  SWCGP. 


2.1.7  Stopping  Criteria  for  Greedy  Algorithms 


Several  criteria  can  be  used  as  halting  rules  for  the  greedy  algorithms.  The  first  option  is  to 
stop  the  greedy  algorithms  when  the  U-norrn  of  the  residual  r^  falls  below  a  preset  threshold. 
This  criterion  aims  to  achieve  a  certain  bound  on  the  reconstruction  error. 

The  second  criterion  is  based  on  the  sparsity  level  of  the  solution  vector  xU,  i.e.,  the  number 
of  nonzero  elements  in  the  solution  vector.  This  stopping  criterion  is  normally  used  when  the 
sparsity  level  of  the  actual  solution  is  a  priori  known  or  a  fixed  number  of  elements  is  desired 
to  approximate  the  solution. 

The  greedy  algorithms  can  also  be  halted  when  the  change  in  the  /2-norm  of  residual  or 
the  maximum  correlation  between  the  residual  and  the  atoms  fall  below  some  threshold. 
These  two  criteria  ensure  the  greedy  algorithms  achieve  a  certain  bound  on  the  reconstruction 
error.  Depending  on  the  applications,  one  or  more  criteria  can  be  applied  to  the  greedy 
algorithms. 
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Algorithm  8  IRLS  Algorithm 

1:  procedure  INPUT:  y,  K.  OUTPUT:  rW,  xU 
2:  Initialisations:  =  <jUy,  =  L 

3:  for  iteration  i  =  1]  i  \=  i  +  1  until  stopping  criterion  is  met  do 

4:  Compute  Weights: 

5:  wf  =  ((|x'l_11|)2  +  tf[<-l])P/2-l  (j  =  1,  .  .  .  ,  N) 

6:  Solution  for  /2-norrn  Objective: 

7:  xW  =  Q[*]$T($Q[il$r)-1y 

8:  where  QM  is  the  diagonal  matrix  with  entries  of  1/w^ 

9:  Updating  the  Regularisation  Parameter: 

10:  if  (||x^  || 2  —  ||x)*— ^ 1 1 2 )  <  V i?M/100  then 

11:  ^[m]  =  ^W/10 

12:  else 

13: 

14:  end  if 

15:  end  for 

16:  end  procedure 


2.2  /i-Minimisation  Algorithms 


The  L -minimisation  approach  with  the  formulation  defined  in  (5)  provides  a  powerful  frame¬ 
work  for  sparse  signal  recovery  with  strong  theoretical  performance  guarantees.  If  the  sensing 
matrix  $  satisfies  a  certain  restricted  isometry  property  (RIP),  a  stable  solution  of  the  sparse 
signal  recovery  problem  is  guaranteed  to  be  obtained  through  the  L -minimisation.  For  the 
case  of  noise-free  signal  recovery  with  £>( y)  =  {x  :  3>x  =  y } ,  the  L -minimisation  can  recover 
any  A-sparse  signal  x  exactly  from  as  few  as  0(I\  log(A/A))  measurements  if  the  sensing 
matrix  $  satisfies  the  RIP  of  order  2 A  with  §2k  <  \/2  ~  1  [17,  18].  For  the  case  of  noisy 
signal  recovery  with  B( y)  =  {x  :  ||y  —  <&x||2  <  e},  the  accuracy  of  the  solution  x  to  (5)  is 
bounded  by 

x  -  x|| 2  <  Co  — — 7=r-^-  +  Cie, 


where 


c  =  pi  ~  (1  ~  V^) 

1  —  (1  +  y/2)  621c 


r  =  .  yi  +  52k 
1  1-(1  +  V2  )52k 


(30) 


if  the  sensing  matrix  $  satisfies  the  RIP  of  order  2 A  with  82K  <  y/2-1  [18].  Here  x^-  denotes 
the  best  A'-sparse  approximation  of  x. 


Along  with  provable  performance  guarantees,  the  l\ -minimisation  formulation  (5)  is  a  convex 
optimisation  problem  which  generally  can  be  solved  effectively  via  any  general-purpose  convex 
optimisation  techniques  [48,  49].  Specifically,  the  Zi -minimisation  can  be  posed  as  a  linear 
program  for  B( y)  =  {x  :  <]>x  =  y},  or  it  can  be  considered  as  a  convex  program  with  a  conic 
constraint  for  B( y)  =  {x  :  ||y  —  ^x|| 2  <  e}. 

In  addition  to  the  L -minimisation  formulation,  there  exists  the  unconstrained  l\ -minimisation 
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formulation  which  has  also  received  great  attention  in  the  literature  (see  [20,  50-53]): 


x  =  arg  min 

xeC" 


(31) 


which  is  commonly  referred  to  as  the  ^-penalty  formulation.  This  ^-penalty  formulation  is 
the  core  of  the  well-known  basis  pursuit  denoising  (BPDN)  algorithm  proposed  in  [53].  In  fact, 
the  use  of  Zi-penalty  has  a  long  history  outlined  in  [20].  Here  the  parameter  A  controls  the 
tradeoff  between  the  approximation  error  and  the  sparsity  of  the  approximation  vector.  In  fact, 
the  ^-penalisation  formulation  (31)  is  equivalent  to  the  l\ -minimisation  formulation  (5)  with 
B(y)  =  {x  :  1 1 y  —  <&x||2  <  e}  for  some  particular  value  of  A  [54],  However,  the  value  of  A,  which 
yields  the  equivalence  between  these  formulations,  is  generally  unknown.  Several  methods  for 
determining  A  are  available  in  [53,  55,  56].  Moreover,  the  L-penalty  optimisation  (31)  is  also 
equivalent  to 

x  =  argmin  ||y  —  <l>x|||  subject  to  j|x||i  <  t,  (32) 

xeC" 

for  appropriate  values  of  A.  Here,  t  is  a  positive  parameter  which  is  inversely  related  to  A. 
In  contrast  to  the  l\ -minimisation  which  is  a  quadratically  constrained  linear  program,  the 
optimisation  in  (32)  is  a  quadratic  program.  This  formulation  in  fact  is  used  in  the  least 
absolute  shrinkage  and  selection  operator  (LASSO)  approach  [52]. 


2.3  /^-Minimisation  Algorithms  with  p  <  1 

The  ^-minimisation  approach  aims  to  minimise  the  /p-norm  (p  <  1)  of  the  estimate  x, 

x  =  argmin  ||x||p  subject  to  x  £  £>( y).  (33) 

xeC" 

In  fact,  the  works  in  [33-35]  focused  on  the  case  of  noise-free  measurements  y  =  3>x,  i.e., 
B( y)  =  {x  :  <&x  =  y}.  It  is  empirically  demonstrated  in  [33]  that  exact  recovery  of  sparse 
signals  can  be  achieved  with  substantially  fewer  measurements  by  replacing  the  Zi-norm  by  the 
Zp-norm  with  p  <  1.  In  addition,  the  theoretical  RIP  condition  of  $  for  the  lp- minimisation 
to  produce  an  exact  reconstruction  of  x  is  [33] 


<W  +  b8(a+ i)A'  <6-1  (34) 

where 

6  >  1,  a  =  ?/A2-p)_  with  p  G  (0, 1].  (35) 

In  fact,  this  is  the  generalisation  of  the  RIP  condition  for  p  =  1  and  6  =  3  presented  in  [17]. 
More  importantly,  this  result  implies  that  a  weaker  condition  for  exact  recovery  is  obtained 
for  smaller  p  [33].  Moreover,  for  the  case  of  random  Gaussian  sensing  matrix  <I>,  the  lv- 
minimisation  with  p  £  (0, 1]  can  recover  exactly  any  A'-sparse  signal  x  with  probability 
exceeding  1  —  l/  ( ^ )  given  that  [34] 

M  >  Ci  (p)  K  +  p  C2  {p)  K  log {N/K) ,  (36) 

where  the  constant  C\  and  C2  are  determined  explicitly  and  are  bounded  in  p. 
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The  Zp-minimisation  with  p  <  1  is  a  non-convex  problem  which  is  intractable  as  a  direct 
problem,  as  described  in  the  literature  [34].  However,  when  it  is  cast  as  a  re- weighted  least 
squares  (IRLS)  problem,  it  can  be  solved  via  an  iterative  approach  [34,  35].  The  main  idea 
behind  IRLS  is  that  it  replaces  the  lp  objective  function  in  (33)  by  a  weighted  /2-norm: 

N 

xW=argminV''  u>^x2  subject  to  <hx  =  y  (37) 

xeCN  ' 


where  the  weights  (j  =  1, . . . ,  N)  are  computed  from  x)*  4  obtained  from  the  previous 
iteration.  The  I2  objective  in  (37)  becomes  the  first-order  approximation  of  the  lp  objective 
in  (33)  when  =  |xy  ^|p_2.  However,  a  regularisation  parameter  is  introduced  to  improve 
the  estimation  performance  [34,  35]: 


(i*jri]i)2+^ 


p/2-1 


(38) 


The  full  description  of  the  algorithm  is  given  in  Algorithm  8.  Note  that  this  ^-regularised 
IRLS  algorithm  becomes  the  FOCUSS  algorithm  in  [57]  if  the  fl-regularisation  strategy  is 
removed.  However,  numerical  results  in  [35]  show  that  the  performance  of  the  -d-regularised 
IRLS  algorithm  is  significantly  superior  over  the  performance  of  the  FOCUSS  algorithm. 


3  Comparative  Performance  Study 

In  this  section,  we  focus  our  attention  on  the  greedy  pursuit  family  and  present  a  performance 
comparison  of  different  greedy  pursuit  algorithms  in  which  the  dictionary  is  built  from  tilted- 
wire  atoms — the  exact  solutions  of  the  tilted-wire  scatterer  model. 


3.1  Signal  Analysis  of  Rotating  Blades 

The  scattered  signal  s(t )  in  the  time  domain  t  received  from  a  rotating  blade  can  be  decom¬ 
posed  as  the  weighted  sum  of  I\  tilted-wire  components  [13,  14], 

K 

s(t)  =  J2Pkg(t,Ok),  (39) 

k= 1 


where 


5(t;i?fc)  =  g{t]rk,^k,Lklak) 

=  A  sine  |  sin(Ot  +  ij)k  +  ak)  1  ex.p{ibrk  sin(Qkt  +  ipk)}  (40) 

is  the  tilted-wire  scatterer  model  of  the  received  signal  from  a  tilted  straight  wire  of  length 
Lk  with  radial  distance  rk  to  its  mid-point  from  the  origin  as  depicted  in  Figure  1,  while  pk 
is  its  (generally  complex- valued)  coefficient.  The  wire  rotates  around  the  origin  with  angular 
velocity  H  from  initial  angle  -*/> k ,  and  is  tilted  at  a  fixed  angle  ak.  In  (40),  A  is  the  normalisation 
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Figure  1:  Tilted  wire  geometry  (radar  is  located  in  the  positive-y  direction.) 


constant,  and  the  parameter  b  is  defined  by  b  =  47r/A  with  A  being  the  wavelength  of  the  radar 
signal. 

In  the  context  of  compressed  sensing,  we  refer  to  g{t\  $&)  as  a  ‘tilted-wire  atom’.  Note 
that  sinc()  is  the  unnormalized  sine  function,  and  i  =  \J—  1.  Writing  (39)  in  vector  form 
yields 

K 

s  =  ^2PkSk,  (41) 

k= 1 

where 

s  =  [s(U),  •  •  • ,  s(tN)]T ,  and  gfc  =  [g(h,  &k), . . . ,  g(tN ,  tfk)]T  (42) 

are  vectors  of  N  discrete-time  samples  of  s(t)  and  g(t,  $&)  respectively.  Since  an  object  can 
be  represented  as  a  small  number  of  scattering  elements,  the  problem  of  radar  imaging  for 
rotating  blades  can  be  cast  into  the  sparse  signal  representation  problem  over  an  over  complete 
library  of  tilted-wire  atoms.  Specifically,  the  aim  is  to  find  a  sparse  solution  of  the  coefficient 
vector 

P  =  \pii  ■  ■  ■  1  Pm]  (43) 

from 

s  =  Gp,  (44) 


where  the  dictionary 


G  =  [gl,...,gAf] 


(45) 


consists  of  M  tilted-wire  atoms  spanning  over  a  discrete  (generally  4-dimensional)  grid  of 
parameters  Note  that,  since  the  scattering  characteristics  from  the  approaching  edges 
(i.e.,  resulting  positive  Doppler  frequency)  are  different  to  those  from  the  receding  edges,  the 
positive  and  negative  Doppler  flashes  are  processed  separately  [14],  This  can  be  done  by 
including  both  positive  and  negative  radius  values  [14] . 
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3.2  Results  With  Simulated  Data 


Consider  a  simulated  helicopter  rotor  system  with  three  blades  rotating  at  40  rad/s,  with  each 
blade  consisting  of  three  wires  with  parameters  given  in  Table  1.  The  parameter  grids  used 
for  constructing  the  function  dictionary  G  are  1  <  L  <  5  nr  in  steps  of  0.5  m,  —  4  <  r  <  4  m 
in  steps  of  0.5  nr,  —1°  <  a  <  1°  in  steps  of  0.5°,  and  ip  spanning  over  a  interval  of  ±5°  in 
steps  of  0.5°  around  —12.5°,  —132.5°  and  107.5°,  which  are  the  angular  positions  ip  of  the 
blades  when  they  are  orthogonal  to  the  radar.  The  wire  parameters  have  been  selected  such 
that  they  are  relatively  very  close  in  the  physical  parameter  space,  resulting  in  a  significant 
level  of  correlation  among  ‘nearby’  atoms  in  the  dictionary. 

Table  1:  Tilted-wire  parameters  used  for  the  target  model. 


Blade 

Wire 

Weighting  Coeff. 

ip  (deg) 

a  (deg) 

r  (m) 

L  (m) 

1 

4.5 

-12 

-0.5 

2.5 

3 

1 

2 

5 

-12.5 

0 

3 

4 

3 

4 

-13 

0 

1.5 

3 

1 

4.5 

-132 

-0.5 

2.5 

3 

2 

2 

5 

-132.5 

0 

3 

4 

3 

4 

-133 

0 

1.5 

3 

1 

4.5 

108 

-0.5 

2.5 

3 

3 

2 

5 

107.5 

0 

3 

4 

3 

4 

107 

0 

1.5 

3 

5 

4.9 
4.8 
4.7 
4.6 
4.5 
4.4 
4.3 
4.2 
4.1 
4 

-6  -4  -2  0  2  4  6  -0.2  0  0.2 

Cross-range  (m)  Cross-range  (m) 

(a)  Orientation  of  the  three  wires  (b)  Exploded  view  of  a  single  wire 

Figure  2:  Scatter  plots  of  the  target  model  showing  the  orientation  and  length  of  the  3  tilted 
wires.  Color  bar  indicates  the  amplitude  of  the  wire  coefficients. 
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Figure  2  is  a  plot  of  the  3-bladed  rotor  showing  the  orientation  and  length  of  each  wire.  The 
synthetic  data  for  the  received  signal,  obtained  by  a  pulse-Doppler  radar  operating  at  the 
transmitted  frequency  of  9.5  GHz  and  the  pulse  repetition  frequency  (PRF)  of  66  kHz,  is 
generated  using  (39)  for  one  cycle  of  rotor  rotation.  The  signal-to-noise  ratio  (SNR)  is  set  to 
20  dB.  The  time-domain  and  spectrogram  plots  of  the  synthetic  received  signal  are  given  in 
Figure  3. 


N 

i  10 

>, 

c  0 
0 

0-10 
I— 

LL 

50  100  150 

Time  (ms) 

Figure  3:  Time-domain  and  spectrogram  plots  of  the  synthetic  signal. 

Figure  4  shows  the  comparison  in  the  time  domain  of  the  reconstructed  signals  obtained 
by  MP,  OMP,  CGP,  gOMP,  CoSaMP  and  ROMP  for  a  single  Monte  Carlo  (MC)  run.  In 
the  frequency  domain,  the  techniques  appear  almost  the  same,  with  no  physically  significant 
differences.  Figure  5  provides  the  corresponding  scatter  plots  of  all  tilted-wire  atoms  found  in 
the  signal  representation.  Taking  into  account  the  relative  magnitudes  of  representing  atoms, 
again  the  techniques  produce  very  similar  results.  Here,  the  results  in  Figures  4-5  are  obtained 
using  12  atoms  per  blade  flash.  These  results  indicate,  at  least  qualitatively,  that  even  when 
the  scatterers  are  very  close  together  in  the  spatial  parameter  space,  the  MP,  OMP,  CGP, 
and  CoSaMP  techniques  still  perform  quite  satisfactorily  at  least  in  term  of  representing  the 
original  signal  in  time  domain. 

Figures  6-8  show  the  root-mean-squared-error  (RMSE)  of  the  reconstruction  signal  and  the 
averaged  running  time  versus  the  sparsity  level  in  the  solution  vector  for  MP,  OMP,  CGP, 
gOMP,  CoSaMP  and  ROMP.  The  RMSE  of  the  reconstructed  signal  is  computed  over  short 
intervals  with  51  samples  around  the  six  main  blade  flashes  and  is  normalised  by  the  U-norm 
of  the  original  signal  in  the  corresponding  intervals.  The  sparsity  level  in  the  solution  vector 
means  the  number  of  collected  atoms  for  each  blade  flash.  Here,  the  results  in  Figures  6-8 
were  obtained  using  50  MC  runs. 

A  general  indication  observable  from  these  results  is  as  follows:  OMP  is  the  most  accurate 
(smallest  reconstruction  error  in  a  least-squares  sense)  but  most  computationally  expensive. 
CGP  is  very  similar  to  OMP  in  terms  of  accuracy  performance,  especially  for  high  sparsity 
levels,  and  with  approximately  one  third  the  computational  cost  of  OMP.  MP  is  simple,  as 
fast  as  CGP,  but  with  an  error  of  about  70  or  80  percent  larger  than  that  of  OMP.  The  other 
techniques  are  as  fast  as  CGP  but  with  excessive  errors,  at  least  for  this  particular  simulated 
target  model. 
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Figure  4-  Comparison  in  time  domain  of  the  reconstructed  signals  (red  trace)  and  the  original 
signal  (blue  trace)  around  the  first  blade  flash  of  the  simulated  data,  for  each  of  the 
greedy  algorithms. 
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Figure  5:  Scatter  plots  of  the  extracted  tilted-wire  atoms  representing  the  simulated  signal, 
with  each  of  the  greedy  algorithms.  Color  bar  indicates  the  atom  amplitude. 
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Figure  6:  Normalised  reconstruction  error  and 
using  MP,  OMP  and  CGP. 
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Figure  7:  Normalised  reconstruction  error  and  running  time  after  processing  synthetic  data 
using  gOMP. 
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Figure  8:  Normalised  reconstruction  error  and  running  time  after  processing  synthetic  data 
using  ROMP,  CoSaMP  and  OMP. 
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3.3  Results  With  Real  Data 

This  section  uses  real  data  from  a  3-bladed  Squirrel  helicopter  collected  with  an  experimental 
radar  operating  at  the  frequency  of  9.5  GHz  and  the  PRF  of  66  kHz.  The  raw  signal  is  high- 
pass  filtered  to  attenuate  low-frequency  components.  The  rotation  speed  and  the  number  of 
rotor  blades  are  estimated  based  on  [58] .  A  rough  estimate  of  the  initial  angles  of  the  blades 
are  also  obtained  based  on  the  location  of  the  main  blade  flashes  in  the  time-domain  signal. 
Figure  9  displays  the  time-domain  and  spectrogram  plots  of  the  preprocessed  signal.  Though 
the  pre-processed  signal  contains  all  return  components  from  the  helicopter  which  include  the 
main  rotor,  tail  rotor,  rotor  hub,  and  the  fuselage  (or  aircraft  body),  only  the  main  rotor 
blades  will  be  analysed  in  this  study. 

The  function  dictionary  used  here  is  constructed  from  a  parameter  grid  for  U  spanning  over 
interval  of  ±5°  in  steps  of  0.5°  around  these  initial  angle  estimates.  The  same  parameter 
grids  for  L,  r  and  a  are  used  as  in  Section  3.2  for  the  simulated  data  example.  Here,  a 
coherent  processing  interval  of  150  ms  is  used  corresponding  to  approximately  one  cycle  of 
rotor  revolution. 

Figures  10  and  11  show  the  time-domain  and  frequency-domain  plots  of  the  reconstructed 
signals  obtained  by  MP,  OMP,  CGP,  gOMP,  CoSaMP  and  ROMP  as  compared  to  the  original 
real  signal.  The  agreement  is  reasonable  though  not  as  close  as  in  the  simulated  example, 
which  can  be  attributed  to  various  real  effects:  a  real  rotor  blade  in  flight  deviates  signifantly 
from  the  straight  and  static  shape  -  it  flaps  and  bends  in  a  rather  random  and  chaotic  manner, 
as  can  be  verified  elsewhere  [59].  Another  interesting  observation  is  that  the  components  in 
lower  Doppler  frequency  regions  near  the  rotor  hub  were  picked  up  by  MP,  OMP,  and  CGP, 
but  not  by  the  other  techniques.  Note  that  MP,  OMP,  and  CGP  pick  up  only  one  atom  in  the 
Identify  step  of  each  iteration,  whereas  the  other  techniques  pick  up  multiple  atoms. 

Figure  12  provides  the  corresponding  scatter  plots  showing  the  wire  model  parameters  for 
all  the  collected  atoms.  Again,  the  results  in  Figures  10-12  are  obtained  using  12  atoms  per 
flash.  The  relative  accuracies  of  these  techniques  can  be  qualitatively  assessed  through  the 
locations  of  the  wire  elements:  MP,  OMP,  and  CGP  produce  reasonable  wire-frame  images, 
while  CoSaMP  is  the  worst. 

The  quantitative  results  in  Figures  13-15  confirm  the  above  qualitative  assessment.  As  before, 
these  plots  show  the  normalised  /2-norm  of  the  reconstruction  error  and  the  running  time 
versus  the  sparsity  level  in  the  solution  vector  for  MP,  OMP,  CGP,  gOMP,  CoSaMP  and 
ROMP.  The  U-norm  of  the  reconstruction  error  is  computed  over  short  intervals  with  51 
samples  around  the  six  main  blade  flashes  and  is  normalised  by  the  ^-norm  of  the  original 
signal  in  the  corresponding  intervals.  The  sparsity  level  in  the  solution  vector  means  the 
number  of  collected  atoms  for  each  flash. 
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Figure  9:  Time-domain  and  spectrogram  (frequency-domain)  plots  of  the  Squirrel  helicopter 
data. 
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Figure  10:  Comparison  in  the  time  domain  of  the  reconstructed  signals  ( red  line)  and  the 
original  signal  (blue  line)  around  the  first  blade  flash  in  real  data,  for  each  of  the 
greedy  algorithms. 
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Figure  11:  Spectrogram  plots  of  the  reconstructed  signals  after  processing  real  data  with  each 
of  the  greedy  algorithms. 
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Figure  12:  Scatter  plots  showing  the  tilted- wire  atom  parameters  after  processing  real  data 
with  each  of  the  greedy  algorithms.  Color  bar  indicates  the  amplitude  of  the  atom 
coefficients. 
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Figure  13:  Normalised  reconstruction  error 
MP,  OMP  and  CGP. 
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Figure  If:  Normalised  reconstruction  error  and 
gOMP. 
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Figure  15:  Normalised  reconstruction  error  and  running  time  after  processing  real  data  using 
ROMP ,  CoSaMP  and  OMP. 
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3.4  Further  Discussion 

The  report  documents  a  comprehensive  survey  of  the  currently  known  sparsity-based  signal 
processing  techniques.  From  the  results  for  both  simulated  and  real  data,  we  observe  that 
OMP,  MP  and  CGP  exhibit  similar  performance  where  OMP  is  slightly  superior  than  MP 
and  CGP.  However,  OMP  takes  much  longer  time  to  run  compared  to  that  of  MP  and  CGP. 
Since  CGP  only  performs  a  single  gradient  update  at  each  iteration,  CGP  requires  a  sufficient 
number  of  iterations  to  converge  to  the  ‘correct’  least-squares  solution.  For  that  reason,  CGP 
performs  worse  than  OMP  for  low  sparsity  levels  of  the  solution  vector.  On  the  other  hand, 
for  sufficiently  large  sparsity  in  the  solution  vector,  CGP  performs  very  similar  to  OMP  while 
it  is  more  efficient  in  terms  of  computation. 

The  results  also  suggest  that,  by  allowing  the  selection  of  multiple  atoms  in  each  iteration, 
gOMP  gains  significant  reduction  in  runtime  compared  to  OMP.  However,  this  computational 
advantage  comes  at  the  expense  of  performance  degradation  as  expected.  A  similar  trend 
was  observed  for  stOMP,  SWOMP  and  SWCGP.  The  performance  of  stOMP,  SWOMP  and 
SWCGP,  however,  is  not  shown  here  as  they  are  similar,  at  least  for  these  examples.  In 
addition,  although  providing  very  strong  theoretical  guarantees,  it  is  observed  that  ROMP, 
SP  and  CoSaMP  do  not  perform  well  compared  to  OMP.  Note  that,  due  to  the  similarity  with 
SP  and  CoSaMP,  only  the  performance  of  CoSaMP  is  included  in  this  report.  The  common 
feature  of  gOMP,  SWOMP,  SWCGP,  ROMP,  SP  and  CoSaMP  is  that  they  select  multiple 
new  atoms  at  a  time.  However,  for  the  particular  problem  of  radar  imaging  of  rotating 
blades,  selecting  multiple  atoms  is  not  desirable  because  the  function  dictionary  G  may  be 
constructed  from  fairly  dense  grids  of  L ,  r,  U  and  a.  As  the  result,  along  with  the  correct 
atom,  these  algorithms  may  select  other  atoms  which  are  closely  located  with  the  correct  one, 
disrupting  the  accuracy  of  the  representation. 

Figures  16-21  showing  the  evolution  of  collected  atoms  with  iteration  index  provide  further 
insight  into  the  behaviour  of  the  techniques.  Here,  the  algorithms  for  MP,  OMP,  CGP,  gOMP 
and  ROMP  terminate  when  10  atoms  are  selected  while  CoSaMP  terminates  when  the  I2- 
norm  of  the  signal  reconstruction  error  dips  below  a  certain  threshold.  Note  that  the  stopping 
criterion  based  on  the  sparsity  level  in  the  solution  vector  is  not  appropriate  for  CoSaMP 
because  it  maintains  a  fixed  sparsity  level  in  the  solution  vector  over  iterations  by  removing 
insignificant  atoms  via  the  pruning  step.  In  Figures  16-21,  the  magnitude  of  coefficient  is 
color  coded  and  the  coordinate  system  has  been  rotated  so  that  the  blade  is  oriented  upward. 
As  shown  in  Figures  16-18,  for  the  greedy  algorithms  (i.e. ,  MP,  OMP  and  CGP)  which  select 
a  single  atoms  at  a  time,  the  coefficients  of  the  collected  atoms  are  very  stable  over  iterations. 
On  the  other  hand,  gOMP,  ROMP  and  CoSaMP  select  a  group  of  atoms  which  are  closely 
located  as  shown  in  Figures  19-21.  As  a  result,  along  with  the  ‘most  matched’  atom,  these 
algorithms  tend  to  also  select  other  incorrect  atoms  ‘nearby’.  This  is  a  plausible  explanation 
for  the  performance  degradation  of  gOMP,  ROMP  and  CoSaMP  compared  to  OMP. 

We  have  also  carried  out  other  simulation  runs  where  the  true  scatterers  are  more  separated, 
the  performances  of  the  techniques  become  more  comparable;  however,  this  does  deviate  from 
the  current  interest  and  will  be  further  investigated  in  future  studies. 
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Figure  16:  Scatter  plots  showing  the  evolution  of  atom  parameters  extracted  using  MP. 
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Figure  17:  Scatter  plots  showing  the  evolution  of  atom  parameters  extracted  using  OMP. 
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Figure  18:  Scatter  plots  showing  the  evolution  of  atom  parameters  extracted  using  CGP. 
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Figure  19:  Scatter  plots  showing  the  evolution  of  atom  parameters  extracted  using  gOMP. 
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Figure  20:  Scatter  plots  showing  the  evolution  of  atom  parameters  extracted  using  ROMP. 
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4  Conclusion 


This  report  is  part  of  a  study  on  techniques  for  the  analysis  of  radar  backscatter  from  fast 
rotating  blade- like  objects  using  sparsity-based  approaches  in  which  the  function  dictionary  is 
based  on  the  exact  solution  of  the  tilted-wire  scatterer  model  in  two  spatial  dimensions.  Three 
classes  of  sparse  recovery  algorithms  were  reviewed,  based  on  the  greedy,  L -minimisation  and 
Zp-minimisation  (p  <  1)  methods.  However,  most  of  the  report  focuses  on  the  performance 
and  computational  cost  of  the  greedy  algorithms. 

Using  both  simulated  and  real  radar  data,  this  preliminary  study  found  that  those  greedy 
algorithms  that  select  multiple  atoms  at  each  iteration  perform  poorly  in  terms  of  high  error 
in  the  signal  reconstruction,  although  the  computational  cost  is  generally  many  times  lower; 
these  include  the  techniques  of  gOMP,  SWOMP,  SWCGP,  ROMP,  SP  and  CoSaMP.  On  the 
contrary,  the  techniques  of  MP,  OMP,  and  CGP  which  select  only  one  atom  at  each  iteration 
perform  better  in  terms  of  representation  error,  although  with  a  higher  computational  cost. 
In  particular,  the  CGP  is  computationally  much  more  efficient  than  the  OMP  and  yet  exhibits 
a  very  similar  performance  to  the  OMP.  Thus,  in  terms  of  performance-versus-computational 
cost  trade-off,  the  CGP  is  the  algorithm  of  choice  for  sparse  analysis  of  narrowband  radar 
backscatter  from  helicopter  blades.  This  conclusion  is  valid  at  least  for  the  cases  where  high 
correlation  exists  among  the  atoms  of  the  underlying  dictionary,  which  arises  from  closely 
spaced  scatterers  in  the  spatial  domain. 

Further  investigation  is  currently  underway  and  will  be  reported  in  future  publications. 
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